#ifndef MATRIX4X4_H
#define MATRIX4X4_H

typedef struct matrix4x4
{
  double data[ 4 ][ 4 ];
} matrix4x4_t;


matrix4x4_t *
matrix4x4_new( void );

void
matrix4x4_free( matrix4x4_t *m );

char *
matrix4x4_str( const matrix4x4_t m );

static inline void
matrix4x4_set_identity( matrix4x4_t *m )
{
  unsigned int i, j;
  for ( i=0; i < 4; i++ )
    for ( j=0; j < 4; j++ )
      m->data[ i ][ j ] = ( i == j );
}

static inline matrix4x4_t *
matrix4x4_new_identity( void )
{
  matrix4x4_t *m = matrix4x4_new();
  matrix4x4_set_identity( m );
  return m;
}

static inline matrix4x4_t *
matrix4x4_multiply( const matrix4x4_t a, const matrix4x4_t b )
{
  matrix4x4_t *r = matrix4x4_new();
  unsigned int i, j, k;

  for ( i=0; i < 4; i++ )
    for ( j=0; j < 4; j++ )
      for ( k=0; k < 4; k++ )
	r->data[ i ][ j ] += a.data[ i ][ k ] * b.data[ k ][ j ];

  return r;
}

static inline matrix4x4_t *
matrix4x4_transpose( const matrix4x4_t m )
{
  matrix4x4_t *r = matrix4x4_new();
  unsigned int i, j;

  for ( i=0; i < 4; i++ )
    for ( j=0; j < 4; j++ )
      r->data[ i ][ j ] = m.data[ j ][ i ];

  return r;
}

static inline matrix4x4_t *
matrix4x4_copy( const matrix4x4_t m )
{
  matrix4x4_t *r = matrix4x4_new();
  unsigned int i, j;

  for ( i=0; i < 4; i++ )
    for ( j=0; j < 4; j++ )
      r->data[ i ][ j ] = m.data[ i ][ j ];

  return r;
}


#endif /* MATRIX4X4_H */
